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We report on the stress response of gyroidal and lamellar amphiphilic mesophases to steady shear 
simulated using a bottom-up lattice-Boltzmann model for amphiphilic fluids and sliding periodic 
(Lees-Edwards) boundary conditions. We study the gyroid per se (above the sponge-gyroid transi- 
tion, of high crystallinity) and the molten gyroid (within such a transition, of shorter-range order) . 
We find that both mesophases exhibit shear-thinning, more pronounced and at lower strain rates 
for the molten gyroid. At late times after the onset of shear, the skeleton of the crystalline gyroid 
becomes a structure of interconnected irregular tubes and toroidal rings, mostly oriented along the 
velocity ramp imposed by the shear, in contradistinction with free-energy Langevin-diffusion studies 
which yield a much simpler structure of disentangled tubes. We also compare the shear stress and 
deformation of lamellar mesophases with and without amphiphile when subjected to the same shear 
flow applied normal to the lamellae. We find that the presence of amphiphile allows (a) the shear 
stress at late times to be higher than in the case without amphiphile, and (b) the formation of rich 
patterns on the sheared interface, characterised by alternating regions of high and low curvature. 



I. INTRODUCTION 

The study of the response to shear in amphiphilic 
mesophases has been the subject of attention for numer- 
ical modellers only in recent years. The interest in the 
subject is sustained not only by the wide range of applica- 
tions in materials' science and chemical engineering, but 
also by the need to gain a fundamental understanding of 
the universal laws governing the self-assembly processes 
and competing mechanisms present. 

Hitherto, studies have focused mainly on the struc- 
tural changes induced by steady and oscillatory shear, 
near and far from critical points, in polymer sys- 
tems [1-6]. The morphologies studied include cubic- 
and wormlike- micellar, lamellar and hexagonally-packed- 
tubular mesophases; more complex structures are the so- 
called bicontinuous mesophases, of which those liquid- 
crystalline of cubic symmetry have thus far been consid- 
ered in far less detail. 

The amphiphilic gyroid [14, 15] is a bicontinuous cubic 
liquid crystal consisting of multi- or mono-layer sheets of 
self- assembled amphiphile dividing two regions, each con- 
taining phases which are mutually immisiciblc, e.g. aque- 
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ous and hydrocarbon species. These sheets or labyrinths 
form a triply periodic minimal surface (TPMS) whose 
unit cell is of cubic symmetry, has zero mean curvature, 
no two points on it are connected by a straight segment, 
and has no reflexion symmetries. Their skeletons, i.e. 
the locus bounded by the TPMS, for each immiscible 
phase, form double (inter- weaving) , chirally symmetric 
3-fold coordinated lattices. There are lyotropic [14, 15] 
and thermotropic transitions between the gyroid and the 
microemulsion mesophase, the latter being a bicontinu- 
ous mesophase of short-range order. The morphologies 
in the crossover regions of these transitions show shorter- 
range order than the gyroid's and longer-range order than 
the microemulsion's, reasons for which they are termed 
'molten gyroids.' 

Bicontinuous cubic mesophases of monoglycerides and 
the lipid extract from archaebacterium Sulfolobus solfa- 
taricus have been found at physiological conditions in cell 
organelles and physiological transient processes such as 
membrane budding, cell permeation and the digestion of 
fats [7]. They can also be synthesised for important ap- 
plications in membrane protein crystallisation, controlled 
drug release and biosensors [8, 9]. 

The purpose of this paper is to report on the re- 
sponse to shear of gyroid (G), molten-gyroid (MG) and 
lamellar (L Q ) amphiphilic mesophases simulated using a 
bottom-up kinetic-theoretic model for fluid flow. The 
model is based on a lattice-Boltzmann (LB) method, 
which has proved to be a modelling tool alternative to 
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and more efficient and robust than sofisticated methods 
based on continuum equations. This LB method ad- 
heres to a bottom-up complexity paradigm [13] in the 
sense that it is simple, fully particulate and no hypothe- 
ses of desirable macroscopic behaviour are imposed on 
the microdynamics — yet, we have shown in the past its 
ability to simulate correct segregation kinetics for immis- 
cible fluids [20] and non-equilibrium self-assembly into 
amphiphilic mesophases [14, 15]. Knowing that such a 
simple model is capable of simulating these kinetic pro- 
cesses from a purely bottom-up dynamics, in this pa- 
per we investigate how hydrodynamic interactions cou- 
ple with self-assembly and modify the stability and mor- 
phology of the mesophases. The novelty of our work also 
rests in the model's capability to reproduce morphological 
transitions without having to assume a macroscopic, free- 
energy model, used in other LB methods [3, 10] to com- 
pute the diffusive currents substantiating self-assembly. 

In addition, since our method models amphiphilic 
molecules as point dipoles — the simplest possible partic- 
ulate model for an amphiphile — the rhcological features 
emergent from it are expected to be universal for a broad 
range of amphiphilic systems. Finally, most of the nu- 
merical studies measuring the stress response of complex 
fluids to shear reported in the literature deal with phase- 
segregating fluids on one side [10], and the more compli- 
cated polymeric [11] and glassy systems on the other [12]. 
In this respect, the present article stands somewhere in 
between these two. 

Our paper is structured as follows. In the next section 
we briefly introduce the model and describe the bound- 
ary conditions for the imposition of shear. In Section III 
we report on simulation data and conclude that shear- 
thinning occurs for both G and MG mesophases, lead- 
ing to a transition to a mesophase consisting of tubu- 
lar and ring-like structures as the strain increases. In 
section IV we reveal how the presence of amphiphile in 
lamellar mesophases induces the formation of rich inter- 
facial patterns surviving shear and allows higher values of 
stress than in lamellar mesophases without amphiphile. 
Finally, we provide our conclusions in Section V. 



II. THE MODEL AND THE LEES-EDWARDS 
BOUNDARY CONDITIONS 

We utilised an existing bottom-up lattice-Boltzmann 
(LB) model for amphiphilic fluids [14, 15], extended to 
simulate shear flow by means of Lees-Edwards bound- 
ary conditions [16]. The model is in turn based on an 
extension made to the Shan-Chen bottom-up LB model 
for immiscible fluids to model amphiphilic-fluid flow, and 
employs 25 microscopic velocities, of speeds 0, 1 and \/2, 
in three dimensions (D3Q25 lattice) [17, 18]. The model 
uses a BGK (relaxation-time) approximation to the colli- 
sion term of the Boltzmann equation for fluid transport, 



which allows us to simulate, for large enough lattices [19], 
the Navier-Stokes (NS) momentum-balance equation in 
the bulk of each immiscible fluid species, namely "oil" (or 
"red", V) and "water" (or "blue", 'b'). The model allows 
the simulation of correct phase-segregation kinetics in the 
absence [20] and presence [15] of a third, amphiphilic 
(surfactant-like, "s") dipolar species. The model controls 
the inter-particle forces between r, b, and s species via 
coupling parameters ((/br, 5bs, 5ss), and transients are con- 
trolled via relaxation times for densities (r b ,r r ,r s ) with 
an additional relaxation time for the orientation of the 
amphiphile dipoles (r d ). In addition, the model simulates 
the noncquilibrium self-assembly and relaxation dynam- 
ics of sponge (L 3 ) and gyroid mesophases [14, 15]. The 
gyroids that it simulates show rigidity, arising from their 
crystalline ordering, which decreases as the concentration 
of amphiphile is reduced; indeed, a lyotropic transition 
causes the correlation length to decrease towards that 
of a sponge mesophase through a molten-gyroid state. 
This idea is central to the work we present here: we shall 
see that the mesophase's crystalline ordering enhances its 
stress response; indeed, we find shear-thinning to occur 
at higher strain rates for gyroids than for sponges. 

The Lees-Edwards boundary conditions (LEBC) were 
originally proposed by Lees and Edwards in the context of 
molecular dynamics simulations [16]. They showed that 
these boundary conditions would give rise to a desired lin- 
ear, wedged velocity profile whilst avoiding the trouble- 
some spatial inhomogeneities appearing when solid walls 
are used to induce the shear flow [21]. A particular 
realisation of the LEBC on a cartesian simulation box 
[0, N x ] x [0, N y ] x [0, N z ] is established by letting the pe- 
riodic images, at iV x < x < 2N X and — iV x < x < 0, 
move parallel to unit vectors ±z, respectively, both with 
speed U. The LEBC, in their original, molecular dynam- 
ics form, are expressed as a Galilean transformation on 
the position (x,y,z) and velocity (£ XJ £y , £z) co-ordinates 
of a molecule, as follows 

x' = xmodN x 

y' = ymodNy 

( {z + A z )modiV z , x > N Xl 
z' = I zmodN z ,0<x<N x , (1) 

[ (z- A z )modiV z , x < 0, 

& = 

r & + u , x > n x , 

& = I & ,0<x<N x , (2) 
[ & - U ,x<0, 

where A z = UAt is the image's shift at time At after the 
onset of shear. 

An implementation of the LEBC on our LB dynamics 
(LB-LEBC) differs from that used in molecular dynamics 
(MD-LEBC) in that the shift A z is not in general a mul- 
tiple of the lattice unit, as Wagner and Pagonabarraga 
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have pointed out [21], and hence an interpolation scheme 
is needed. This interpolation scheme streams the am- 
phiphile dipoles d(x) and mass densities n£(x) located 
at position x on the shearing wall, where c& is the rele- 
vant discrete molecular velocity, k = 1, . . . , 25, for each 
(fluid and amphiphilic) species a. 

In our LB-LEBC, while the spatial displacement fol- 
lows Eqs. (1), the velocity shift cannot be enforced by re- 
placing the continuum velocity component £ z in Eqs. (2) 
with the discrete microscopic speeds Cfc-z, since the veloc- 
ities Cfc are constant vectors. Instead, this acceleration is 
enforced on the macroscopic fluid velocity around which 
the local-Maxwellian distributes molecules at equilibrium 
and towards which the BGK scheme relaxes, similarly to 
how immiscibility forces are implemented [15, 20]. This 
procedure guarantees that all accelerations in the fluid 
are ruled by the same BGK process, controllable via the 
shape of the distribution function and the relaxation-time 
parameter, including the acceleration due to the shearing 
walls. 

The MD-LEBC give rise, at steady state (late times), 
to a shear state which is Galilean-invariant, i.e. no partic- 
ular plane in the system is favoured over another. This is 
a sine qua non for any shearing method, and our method 
satisfies it too. As regards the unsteady, transient initial 
states, the MD-LEBC are unphysical since they cannot 
provide the molecular specificity (e.g. wall roughness) 
required in an atomistic approach to boundary effects, 
such as density layering and slip at wall. However, meso- 
scopic methods — LB is one of them — in general only de- 
scribe low wavenumbers and frequencies, which means 
that, with respect to MD, (a) the atomistic detail of the 
shearing walls is largely coarse-grained and (b) the fluid 
structure and dynamics are much less sensitive to the 
atomistic detail of the walls. Since most boundary ef- 
fects present in MD are absent in LB, the fact that the 
LE boundary conditions eliminate them does not pose 
a problem. This should be taken with a caveat: our gy- 
roidal mesophases melt when placed in a solid box, which 
means that the approach to equilibrium is sensitive to 
momentum transfer with the walls, and therefore the LE 
boundary conditions do not mimic shearing a mcsophasc 
in confinement. (To our knowledge, no bottom-up sim- 
ulations have ever reported mesophase self-assembly in 
confinement.) Rather, the LE boundary conditions in a 
LB model mimic shearing with walls which are far enough 
from the locations in the system where observables are 
probed such that microscopic boundary effects are ab- 
sent. 

Our LEBC implementation is embedded within an effi- 
cient parallel LB algorithm [22] which allows us to employ 
large lattices and hence reach the small Knudscn num- 
ber limit where (a) regions away from interfaces satisfy 
the incompressible NS equation in the limit of low Mach 
numbers (Ma) [20] , and (b) observables vary by less than 
10% when the lateral lattice dimension is doubled. We 



have previously found that the lattice size guaranteeing 
condition (b) is 128 3 for the parameters generating the 
mesophases investigated here [14, 15]. 

III. SHEARING GYROIDAL MESOPHASES 

We sheared two gyroidal mesophases differing in the 
amount of amphiphile present and the value of the inter- 
amphiphile interaction coupling parameter. Each of these 
structures was allowed to self-assemble from homoge- 
neous mixtures of oil, water and amphiphile using pe- 
riodic boundary conditions. They have been appropri- 
ately characterised by probing direct and Fourier-space 
late-time snapshots of the density order parameter <fr = 
p ml — p watcI - more precisely, they correspond to gyroid (cf. 
Fig. 5(a)) and molten gyroid mesophases, as previously 
reported by us [14, 15]. 

The common parameters used for both gyroids were 
oil and water densities flatly distributed in the range 
< n<°) b = n(°) r < 0.7, coupling strengths g hl = 0.08, 
9bs = —0.006, relaxation times t° = r r = r s = r d = 1, 
and, for the amphiphile's dipoles, (3 — 10 and do = 1. 

Their differing parameters were surfactant densities, 
flatly distributed in the initial homogeneous mixture, in 
the ranges < n^ s < 0.9 for the gyroid and < n^ s < 
0.6 for the molten gyroid, with coupling strengths g ss = 
—0.0045 for the gyroid and g ss = —0.003 for the molten 
gyroid. These values for the gyroid are 50% higher than 
those for the molten gyroid. 

While the gyroid relaxes to a highly crystalline struc- 
ture [23], the molten gyroid shows both shorter-range 
order and stronger temporal fluctuations than the for- 
mer [15]. In order to obtain a sufficiently relaxed molten 
gyroid as an initial condition for the shear, we took the 
structure as evolved up to time step 32 500; regarding 
the gyroid, the time slice chosen was time step 15 000. 
For practical reasons, instead of letting the molten gy- 
roid self-assemble starting from a homogeneous initial 
mixture, we upscaled a smaller molten gyroid, previously 
self-assembled using the same parameters on a 64 3 lat- 
tice [15], to a 128 3 lattice. Upscaling consisted in replicat- 
ing identical copies of the system: the periodic boundary 
conditions used to generate the 64 3 system (a) guarantee 
that the density field is smooth across the replica bound- 
aries, yet, for this same reason, (b) produce a molten 
gyroid with an additional, undesirable long-wavelcngth 
fluctuation whose periodicity is half the lattice size. The 
amplitude of this undesired long-wavelength fluctuation 
relaxes in time to a vanishingly small value, fact which 
provides us with the 128 3 mesophase we seek. We ob- 
served, however, that this relaxation takes place in less 
than 1 000 time steps [23], i.e. it is a fast transient which, 
therefore, does not affect the shear response at the late 
times that we are interested in. In other words, the late- 
time shear response is insensitive to a small perturbation 



4 



in the initial condition. This allowed us to take the up- 
scaled, unrelaxed structure as the initial condition for the 
molten gyroid. 

It is worth noting that we did not require an elongated 
aspect ratio for the lattices along the direction parallel 
to the translation of the shearing walls since spatial den- 
sity fluctuations were much smaller than the lattice size. 
This is not the case when shearing phase-segregating flu- 
ids without an amphiphilic, growth-arresting species, as 
has been previously reported using LB lattices of up to 
128 : 128 : 512 sizes and aspect ratio [24]. 



A. Stress response and transients 

Shear thinning is said to occur when the shear viscosity 
drops as the strain rate increases. For structured fluids 
such as those we study in this paper, the dynamic shear 
viscosity, 77, is not expected to be a constant of the strain 
rate 7 = \{d x u z + d z u x ) as is true of Newtonian fluids, 
for which, 



only nearest neighbour interactions are being considered, 
the virial term reduces to 



P xz =±2?77, 77^77(7). 



(3) 



Here P xz is one off-diagonal component of the pressure 
(or stress) tensor, and the sign, by convention, indicates 
that the pressure is exerted by the fluid element on the 
surroundings ('+') or from the latter on the former ('—'), 
respectively. We adhere in this paper to the second case. 
In our simulations, we apply the steady shear described 
in Section II, i.e., the shear is generated by the two image 
cells of the LB lattice located along the x-axis moving in 
opposite directions. As a consequence, d x u z becomes the 
only non- vanishing component of the velocity gradient, 
which is also true for the P xz component of the stress 
tensor (and P zx , since the physical requirement that the 
vorticity, W = \{d x u z — d z ii x ), remains upper bounded 
requires the stress tensor to be symmetric). 

As we have likewise done previously while computing 
diagonal components of the pressure tensor [14, 15, 20], 
here we measured P xz from its definition as the sum of a 
kinetic term plus a virial mean-field term accounting for 
interactions and giving rise to non-ideal gas behaviour, 
namely, 



p w = EE^( x )( c *- u ( x ))( c *- u «) 

a k 

\ X>* E [v- Q (x)v> a (x) + r(*w(*) 



+ 



(x-x')(x-x') 



(4) 



where ip has the form ip = 1— exp[— n(x)], which saturates 
at high density values in order to avoid unbounded inter- 
particle forces whilst reproducing a meaningful equation 
of state [15]. Since the interaction matrix {g a a} is sym- 
metric with all diagonal elements identically zero, and 



^ E 9aa E ^W^t* + c k) c kC k ■ 

ol^o. k 



(5) 



In the incompressible, low Mach number limit, our LB 
model reproduces the NS equation away from inter- 
faces [17, 19], which describes a Newtonian fluid with 
a viscosity being a well known function of the relaxation 
time. The presence of an interface, characterised by an 
interfacial tension and a bending rigidity, however, intro- 
duces anisotropics in the fluid's stress tensor which can be 
accounted for by a tensorial effective viscosity. Since the 
interface may move, at a speed growing with the strain 
rate, these anisotropies can become unsteady. Our aim is 
then to measure how this viscosity evolves with the strain 
and the strain rate. 

In order to probe the function 77 = 77(7) for both gy- 
roidal mesophases, we measured P xz for a number of 
different applied shear rates. The chosen values for U 
were such that they remained within the incompressibil- 
ity limit, i.e. small compared to the speed of sound on 
the D3Q25 lattice, c s = 3~ 1//2 « 0.58. Values chosen were 
U = 0.05,0.10,0.15,0.20, corresponding to Mach num- 
bers Ma = U/c s = 0.086, 0.17, 0.26, 0.34, respectively. All 
observables we report in this paper are spatial averages, 
at least oni = const, planes where a simple fluid un- 
der the same shear would show translational symmetry 
for the velocity field, i.e., perpendicular to the velocity 
gradient. Since, for reasons of computational cost, we do 
not perform averages over the seed used to generate the 
pseudo-random initial configuration mimicking a homo- 
geneous ternary mixture, we do not provide error bars 
around averages. 

Figure I shows the profile of the stress, for the sheared 
gyroid, along the applied velocity gradient direction. Sev- 
eral curves therein depict the transport of momentum to- 
wards the core (i.e., the plane x = 64) as the strain grows 
as a function of time. Distinctively, the profiles have spa- 
tial fluctuations, a consequence of the gyroid's convoluted 
structure whose interfacial tension locally modifies the 
viscosity expected for a simple fluid. The u z component 
of the velocity field, shown in Fig. 2 and averaged in the 
same way as stated for (— P xz ) in the caption of Fig. f , is 
however not inhomogeneous but follows a transient sim- 
ilar to that expected for a simple fluid: we observe the 
setting up of a steady, smooth and wedge-shaped pro- 
file, except at the borders. Figure 2 also includes the 
behaviour of the averaged velocity profile for the molten 
gyroid MG at late times, and is seen to match that of the 
gyroid G. 

Remaining with the G mesophase, we show in Fig. 3 the 
temporal evolution of the stress displayed in Fig. 1; the 
values plotted are averages of the latter on the 8 < x < 
N x — 8 = 120 interval, which amounts to averaging over 
the whole lattice except thin slabs adjacent to the bound- 
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aries. In addition to Fig. 1, we include higher and lower 
shear velocities, namely U — 0.05, 0.15, 0.20. Note that 
the time evolution of the averaged stress is a succession of 
peaks and troughs, denoting successive intervals of yield 
and recoil, which is a canonical feature of viscoelastic be- 
haviour. Were the strain rate at which the gyroid deforms 
coincident with the applied shear rate, these curves would 
imply shear thinning. In fact, while the increments in ap- 
plied shear rate between these curves are kept constant, 
the increments in the (absolute) values of the stress at 
late times do not remain so but decrease. In Fig. 4 we 
show the stress averaged over time steps 24 000 to 28 000, 
plotted against the true strain rate, where the latter was 
measured from the linear velocity profile generated at 
At > 9 000 (t > 24 000), as displayed in Fig. 2. Figure 4 
clearly shows shear thinning: the slope, i.e., the effective 
viscosity r] cS = dP xz /dj, decreases with the strain rate. 
Figure 4 also contains the analogous curve for the molten 
gyroid, which shows shear thinning for the latter at lower 
strain rates than those at which the gyroid does, and at 
higher intensity, i.e. 

<9?7 off _ dif s 



dj 



< 



molten 



dj 



ryroid 



< 0. 



(6) 



This is the first indication of shear thinning reported by 
means of a bottom-up kinetic-theoretic model for fluid 
flow. 



B. Morphological transitions 

Figure 5 shows the configuration of the gyroid in the 
40 < y < 52 slab of the 128 : 128 : 128 lattice, before and 
at late times after applying a shear of U = 0.20. The 
volume-rendering graphical representation employed [25] 
makes regions where <f> > 0.37 opaque to the lighting 
rays, assumed to shine normal to the plane of the text 
and inwards; since —0.79 < 4> < 0-79 over the entire 
system, these regions are the high-density locus of one of 
the species (say, oil). Before shear, the structure contains 
highly ordered subvolumes of gyroid symmetry and diago- 
nal length from about 32 to 64 lattice sites, cf. Fig. 5(a). 
This gyroid is hence a collection of subvolumes with a 
regular tubular structure making up two three-fold co- 
ordinated, interweaving chiral lattices of which we depict 
only one. Since the size of the G unit cell is approximately 
5 to 6 lattice units, the depth (y-dimension) of the slabs 
shown in Fig. 5 is of about two gyroid unit cells. As can 
be seen in Fig. 5(a), the interfaces between these gyroid 
subvolumes are defective regions where long-range order 
and symmetry appears to be drastically reduced [14, 15]. 
Two features characterising them are the spatial varia- 
tion in coordination number and chirality, seen by the 
presence of elongated tubules and toroidal rings, cf. fig- 
ure 6. 

At At = 21 000, which is a late time after the onset 
of shear and we take as steady state, the structure has 



lost any resemblance with the initial gyroid, except for 
the persistence of the toroidal rings, see Fig. 5(c), which 
are defects in G. Also, the structure at At = 21 000 is 
essentially the same as that at time step At = 5 000 — it 
is a nonequilibrium steady state for at least the previous 
16 000 time steps, a time longer than that required for the 
initial configuration to self-assemble from a homogeneous 
mixture of oil, water and amphiphile. The structure at 
At = 21 000 consists of an irregular network of mainly 
the same structural elements characterising the defective 
regions before the onset of shear, namely, (a) elongated 
tubules, with a tendency to align along a direction which 
is a linear combination of directions (1,0, 0) and (0, 0, 1), 
and (b) toroidal, ring-like structures. This description is, 
by visual inspection, similar for every subvolume of the 
lattice visualised. 

We also looked into the structure of the sheared molten 
gyroid at late times. In contradistinction to the gyroid's 
state at high strain, showing tubules of shape similar to 
that depicted in Fig. 6 and at an angle with the x = const, 
planes, the highly strained molten gyroid displays tubes 
which are more stretched and aligned along the z direc- 
tion. The toroidal rings, also present for the molten gy- 
roid before shear, represent a much smaller volume frac- 
tion for the sheared molten gyroid than for the sheared 
gyroid. 

Figure 7 shows the summed structure function 
S(k,t), or scattering pattern, of the sheared gy- 
roid mcsophase, showing stages of its plastic deforma- 
tion. Here, S(k, t) is the structure function, computed 
according to [15, 20] 



S(k,t) = ^ 



Ait) 



(7) 



Here, k is the discrete wavevector, V is the lattice volume, 
C is the unit cell volume for the D3Q25 lattice, and (j>' k (t) 
is the Fourier transform of the fluctuations of (f>. S(k, t) 
is the Fourier transform of the autocorrelation function 
for the order parameter, 



C w (r,t) = (0((x,i)Mx + r,t)) 



(8) 



where r is a vector lag and the brackets indicate an av- 
erage over the spatial coordinate x. Figures 7(a), (b) 
and (d) are the xz 'scattering patterns' of the structures 
in Fig. 5, produced by summing up the structure func- 
tion along the x direction. At At = 1000 (not shown), 
the maximum intensity is reduced to 29% of its value 
at At = 0, while there appear horizontal 'smeared out 
filaments' of very weak intensity, intrinsically related to 
the shearing process, as we shall conclude from Fig. 8. 
At At = 5 000 a clear cardioid shape has developed; the 
fact that it persists for the rest of the simulation confirms 
our observation that the system reaches a steady state at 
time step At = 5 000. In addition, there is no trace of 
gyroidal patterns along the x-direction. 
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In order to investigate the origin of the cardioid shape, 
we computed the scattering pattern for a 'synthetic gy- 
roid', 

G(x) = sin qx cos qy + sin qy cos(qz — S(x) ) + 

sin(gz — <5(x)) cos qx . (9) 

where <5(x) = (x — N x /2)S max is a spatially-varying de- 
phase used to obtain a linear strain on the morphology 
(its maximum value, S max , is reached at the lattice bound- 
aries), and q = const, is a wavenumber controlling the 
size of the surface's unit cell. It is known that G(x) = 
for 6 max = is a good approximation to the Schoen "G" 
triply periodic minimal surface of laid cubic symmetry 
referred to as 'the ideal gyroid' hereafter [26]. Figure 8 
shows the scattering patterns for the unstrained morphol- 
ogy and for dephases <5 max = 8, 16. 

Comparing structure function maps in Figs. 7 and 8, at 
the same value of the strain rate, proves useful. For the 
synthetic gyroid, the strain is controlled by the number 
of unit cells that the dephase causes the structure to shift 
at the lattice boundary, following a linear profile as we 
approach the other boundary going through zero strain 
at the lattice core. For our simulated amphiphilic gyroid, 
however, the strain does not follow a linear profile at early 
times; instead, the strain at time t would need to be 
computed from the integral f f x dt'dx <9 x tt z (x, t'), 
where t' is the time parameter. For the purposes of this 
paper, however, such an analysis would be superfluous; 
in fact, Fig. 8 already provides us with enough informa- 
tion to understand the origin of the cardioid shape. For 
all panels, (a), (b) and (c) therein, the position of the 
peaks at k z = (k x /(2n/N) w -14, 15, where N = 128) 
are invariant under the strain (dephase); not so with the 
peaks at k z ^ 0, which shift leftwards. (The shift would 
be rightwards were d x u z < or <5 max .) The shape of 
the maps in Figs. 7(c) and 7(d) is that of a transformed 
scattering pattern shifted leftwards. This transformation 
occurs early, between At = and At = 3 000, and is char- 
acterised by two (strong, S > 700) peaks similar to those 
of the gyroid at k x = 0, and two (weaker, 200 < S < 700) 
peaks at k z = 0. 

IV. A SIMPLER CASE: SHEARING THE 
LAMELLAR MESOPHASE 

In the last section we reported on the gyroid displaying 
higher shear stress than the molten gyroid. Since the 
structural transition between these two mesophases can 
be driven by both the amphiphile density and the inter- 
amphiphilc coupling parameter, as we have reported in 
the past [15], our aim in this section is to elucidate the 
role of the amphiphile density alone on the stress response 
to shear; we choose the lamellar mesophase as the subject 
of study, since this is the mesophase with the simplest 
possible internal interface. 



The initial configuration employed was a cubic 128 3 
lattice with 16 lamellae, stacked perpendicularly to unit 
vector z. The lamellae were of alternating, oil- water com- 
positions, separated by a thin monolayer of amphiphile; 
the thickness of the immiscible and amphiphilic lamel- 
lae were 7 and 1 lattice sites, respectively. We popu- 
lated each lattice site with a value of density kept con- 
stant over the region corresponding to a given species; 
each microscopic velocity is assigned the same fraction 
of this value. We gave amphiphilic regions the densities 
n (°) s = o, 0.80, 0.95, and oil and water regions the densi- 
ties n(°) r = rj(°) b = 0.7. Shear was applied perpendicular 
to the lamellae with the same LEBC employed in the last 
section, with speed U = 0.10. 

Before the onset of shear, the case without amphiphile 
for the lamellar initial condition just described is, a pri- 
ori, a metastable state in our LB model. In fact, the 
structure has a stationary morphology since short-range 
oil-water forces and the absence of fluctuations main- 
tain immiscibility, i.e. a value for the interface steepness, 
|V</f>|; however, a large enough perturbation in 4> may al- 
low a fluctuation in surface tension which drives the entire 
interface to a radically different shape. Another factor 
disrupting this lamellar morphology is shear, which may 
work against the intcrfacial tension by reducing |V</>|; this 
can lead to miscibility (0 = 0) for high enough strain 
rates. Despite these arguments, we observed stability for 
the sheared lamellar mesophase without amphiphile, as 
we report next. 

Figure 9 shows the stress as measured in the same 
fashion performed on the data plotted in Fig. 3, for sev- 
eral amphiphile densities. The behaviour observed is di- 
verse. For zero amphiphile concentration (solid curve), 
the stress reaches a peak at early times before it pro- 
ceeds to a second, lower maximum at late times, going 
through a trough at intermediate times due to the fact 
that | V0| experiences a transient decrease. 

The high-density regions of one of the immiscible 
species (say, oil) is shown in Fig. 10(a) at late times, At = 
8 000; these are representative of the shape of the oil- 
water interface. Away from the boundaries (x = 0, 128), 
there is a large intcrfacial area with zero curvature, where 
we define the curvature as H = d zz x < p(z), x^(z) being 
the curve resulting from projecting the ^ = 0.18 surface 
onto the xz plane. Curiously, we observe three changes 
of curvature as we follow the curve x^z) for y = const., 
namely, H < 0, H > 0, H < 0, H > 0; instead, we would 
have expected the steady, late-time configuration for the 
sheared lamellar mesophase to rather minimise the in- 
terfacial area, leaving only one inflexion point. For fluid 
regimes under conditions of local thermodynamic equi- 
librium, we can think of H 2 as an interfacial free energy 
density associated with the curvature [27] ; in this case, we 
would have expected the steady, late-time configuration 
to also minimise the intcrfacial free energy. 

The stress curve corresponding to n^ s = 0.80 
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(cf. Fig. 9) shows the absence of large troughs, as occurs 
for the n(°) s = case, despite the fact that interfacial 
tension is drastically reduced by the presence of the am- 
phiphilc. In addition, the stress grows at late times to 
higher values than those achieved by the n^ s = case. 
The late-time order-parameter configuration is displayed 
in Fig. 10(b), showing a rich interfacial pattern. Using 
the same arguments as those of the last paragraph, this 
structure could be characterised by a higher curvature en- 
ergy, J d 2 xH 2 , where d 2 x is a measure on the oil- water 
interface, and H is now defined as the inverse radius of 
curvature, parametrised on the arclength, s. Figure 10(b) 
shows similar regions of high curvature at an equal dis- 
tance from the shearing walls, where u z = const., which 
we shall call nodal planes. Also note that the interface, as 
approximately depicted by the boundary of the (f> > 0.22 
volume, joins the lattice boundary at an angle close to 90 
degrees. 

The stress curve for the n^ s = 0.95 case shows a dra- 
matically different situation for the first 5 000 time steps: 
the presence of a trough, deeper than that present for 
the n(°) s = density. After that, there appears a shoot- 
off whereby the stress rapidly grows and equals the late 
time value achieved in the n^ s = 0.80 case, while the 
order-parameter displays a configuration analogous to the 
n (0)s _ o.80 case, cf. Fig. 10(c). By looking at the am- 
phiphile density field, p s (x), for the case n^ s = 0.95, we 
observed that the high curvature regions arise close to 
the boundaries first (At < 1000), and then rapidly move 
away from them as the strain progresses. 

V. CONCLUSIONS 

In this paper we have reported on the shear stress 
response of two gyroidal cubic amphiphilic mesophases 
previously self-assembled using the same bottom-up LB 
model we employ here, namely, the gyroid per se, G, 
which shows high crystallinity at late self-assembly times, 
and the molten gyroid, MG, endowed with shorter-range 
order and located within the sponge-gyroid lyotropic 
structural transition [15]. Shear was imposed via sliding 
periodic (Lees-Edwards) boundary conditions, and we in- 
vestigated the response to several values of the strain rate. 
In addition, in order to investigate the dependence of the 
shear stress on the amphiphile density, we also sheared 
a lamellar mesophase, of much simpler morphology than 
the gyroidal mesophases. 

We found that the gyroidal mesophases exhibit shear 
thinning, more pronounced and at lower strain rates for 
the MG mesophase than the G mesophase. In other 
words, momentum introduced into the system due to 
shear is transported more easily for the mesophase con- 
taining more amphiphile, with longer-range ordering, i.e. 
the effective viscosity is higher for the G mesophase. 

We also found a shear-induced transition from an initial 



gyroidal morphology (G and MG) to a mesophase char- 
acterised by coexisting elongated tubules and toroidal, 
ring-like structures. The features of this mesophase is in 
contrast to those of the mesophase reported by Zvelin- 
dovsky et al. using free-energy Langcvin-diffusion meth- 
ods by shearing a bicontinuous structure reminiscent of 
a molten gyroid [3]. The structure they found is of a 
shorter-range ordering than that of the MG mesophase 
described here, and the high-strain structure consists of 
coexisting lamellae and hexagonally packed tubes elon- 
gated along the direction of the imposed shear veloc- 
ity. Our sheared mesophases also show enlongated tubes 
along this direction, but the structure is far more com- 
plicated than that found by Zvelindovsky et al. in that 
it exhibits remnant toroidal rings and 'hard shoulders' 
reminiscent of gyroidal skeletons; hexagonal packing and 
coexisting lamellae are, on the other hand, absent. 

The shear performs a plastic deformation which ef- 
fectively breaks the links of the gyroidal skeleton; this 
happens as these links interpose an (oil-water) interface 
whose normal, n, is parallel or anti-parallel to the flow 
streamlines, u. In other words, shear effectively applies 
a 'mixing' force which is in competition with the inter- 
particle forces keeping the mesophase in place, namely, 
those controlled by coupling parameters <?br, gbs and g ss . 
Our hypothesis is that adsorbed dipolcs sitting on in- 
terfacial regions at an angle 9 = Z(u, n) other than 
9 = 0, 7r require more work from the shear forces to 
be drawn away from the interface than those regions on 
which the streamline impinges normally, since the mixing 
force goes as cos 9. In particular, since the mixing force 
vanishes for 9 = tt/2, considerably longer interfaces can 
survive the flow — shear induces a preferential direction 
along which the long-range order present before the on- 
set of shear is not reduced. This explains not only the 
formation of the elongated tubules but also their recon- 
nection (increase in coordination number). In fact, the 
toroidal, ring-like structures are not only vestigial gyroid 
defects which have survived the gradient Vu, but are 
also born anew as a result from reconnections. It is rele- 
vant to point out that Padding and Boek, using a coarse- 
grained molecular dynamics model for wormlike micelles, 
reported on the formation of rings when applying steady 
shear to a wormlike micellar mesophase [6] — this is an 
'amphiphile-in-water' binary mesophase, in contrast to 
the 'oil-amphiphilc-water' ternary mesophases we study 
in this paper. 

By applying shear to a lamellar mesophase we found 
that the presence of amphiphile on the oil-water interface 
of the mesophase causes the interface to fold into a wealth 
of structures with a (discrete) translational symmetry 
on planes equidistant to the shearing walls and along 
the direction of the shear velocity. In other words, the 
inter-amphiphile force couples the adsorbed amphiphilic 
dipoles so that the interface locally increases its curvature 
energy density. It is worth investigating whether this 
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local increase is due to the amphiphile being incapable 
of sustaining interfacial regions of low curvature under 
shear, i.e. whether the 'breaking' mechanism induced 
by shear is counteracted by regions of high curvature en- 
ergy density. Regarding the shear stress, our amphiphilc- 
containing lamellae responded with higher stress at late 
times than those without amphiphile. This contrasts with 
the results found for the gyroidal mesophases, and lets us 
conclude that it is the gyroid's cubic morphology that 
allows this structure to be stiffcr. Understanding the be- 
haviour of the lamellar mesophasc under shear requires 
the study of amphiphile self-assembly under shear, in- 
cluding in and out of plane amphiphilic and associated 
Marangoni currents, and their coupling to the imposed 
flow. 
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FIG. 1: Shear stress response of a gyroid mcsophase along the 
direction of the velocity gradient. As initial condition, we have 
taken a gyroid on a N x NyN z = 128 cubic lattice at time step 
t = 15 000 of self-assembly [14, 15]. The Lees-Edwards walls move 
with speed U = 0.10 (Ma = 0.17). For each x coordinate, the 
original field has been averaged on the plane [1, N y ] X [16, N z — 16], 
where the excluded interval on the 2-axis accounts for wrapped- 
round densities. Standard errors of the averages are about 6 X 10~ 8 
throughout, and are not shown. Each line represents the response 
at At time steps after the start of steady shear: At = (dotted 
line), At = 100 (dash-dotted), At = 800 (dashed) and At = 9 000 
(solid), where the last is ca. the time at which the core (i.e., the 
plane x = 64) fully responds. From the figure we can see that 
momentum transfer decreases as it reaches the core from the walls. 
Also, note that the stress inverts its sign at late times adjacent to 
the boundaries, \x — x$\ < 2 (xq = 0, 128). All quantities reported 
are in lattice units. 
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FIG. 2: Spatially averaged velocity component u z for the molten 
gyroid and the gyroid mesophases sheared with U = 0.10, at late 
times and over the x > 64 half of the system. The dashed thin 
and thick curves correspond to the molten gyroid at time steps 
At = 9 000 and 13 000, respectively. The solid thin and thick curves 
correspond to the gyroid at time steps At = 9 000 and 13 000, 
respectively. The average is over the same two-dimensional domain 
as described in Fig. 1, for each x, and its standard error is shown 
as negligible error bars. Note that the velocity shows a maximum 
located from 2 to 4 sites away from the boundary, unlike a simple 
fluid which would display it exactly at the boundary. The value of 
this maximum coincides with the actual velocity at which the BGK 
relaxation process of our LB model is forcing the fluid to move, 
which needs not coincide with the input parameter U = 0.10. Note 
that the inversion in the sign of the stress that we reported in 
Fig. 1 occurs precisely for \x — xo\ < 2, xq = 0, 128 and at (late) 
times close to and after At = 9 000. The behaviour at the other 
boundary region is similar and symmetric to that displayed here. 
All quantities reported are in lattice units. 
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FIG. 3: Temporal evolution of the average shear stress of the 
gyroid for different values of steady shear. The initial condi- 
tion is the same as that mentioned in Fig. 1. The curves, as 
seen, e.g., at At = 4000 from bottom to top, correspond to 
Lees-Edwards walls moving with speeds U = 0.05, 0.10, 0.15, 0.20 
(shear rates S/IO" 3 = 0.39, 0.78, 1.17, 1.56), respectively. The 
dotted curves are the responses after a sudden termination of 
shear; they are also referred to as the system's relaxation func- 
tions for the relevant shear speeds. The average here is in 
the three-dimensional domain [8, 7V X — 8] X [1, N y ] X [16, N z — 16], 
where 7V X = N y = N z = 128 and error bars are the standard error 
of the average. All quantities reported are in lattice units. 
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FIG. 4: Both the gyroid (solid line) and the molten gyroid 
(dashed) mesophases exhibit shear thinning. Shown is the stress 
averaged over the interval 24 000 < t < 28 000. From the figure it 
is clear that the gyroid manifests greater stiffness than the molten 
gyroid and its (effective) viscosity drops for higher strain rates. All 
quantities reported are in lattice units. 
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FIG. 5: High-density locus of one species (say, oil) in the gyroid mesophase, before shearing, (a), and at an early, (b), and late 
time slice, (c), after the onset of shear. The shear speed is U = 0.20. The complementary immiscible fluid (water) fills the 
voids with a similar, inter-weaving structure. The system is on a 128 x 128 x 128 lattice, and all figures show the subvolume 
40 < y < 52 and the reference system in use (the y-axis is perpendicular to the plane of the page). The initial configuration, (a), 
is a gyroid at 15 000 time steps of self-assembly under periodic boundary conditions. These images are volume renderings of the 
density order parameter, 4> = p oA — p wator ; the regions visible to the reader are those for which <f) > 0.36 whilst over the entire 
fluid —0.79 < (f> < 0.79; the regions for which 4> < —0.36 (water, not shown) display a similar structure which is complementary 
(interweaving) to the one shown here. All quantities reported are in lattice units. 




FIG. 6: Schematic representation of the skeleton (locus of 
highest density) of the gyroid mesophase we employ, and two 
of its structural features before and at late times after the on- 
set of steady shear. The thickness provides a sense of perspec- 
tive, and represents how close each segment is to the reader; 
note that the figures on the right are planar. The skeleton 
denoted by 'gyroid' depicts a portion of one of the two chiral 
lattices making up the long-range order regions of the gyroid 
before shear, cf. Fig. 5a — the coordination number is three at 
each node. In the regions of the gyroid containing defects, 
as well as in most of the sheared mesophase at late times, 
the coordination number can be reduced to two, describing 
a 'tubule'. We also show the skeleton of the 'ring' structure 
ubiquitous in the sheared gyroid at late times, also present 
in smaller proportion as a defect in the mesophase before the 
onset of shear. At lower values of density, this ring appears 
to be toroidal. 
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FIG. 7: Projected structure function ('scattering pattern') as a function of the time step for the sheared gyroid, as calculated 
using Eq. (7). Shear velocity is U = 0.20. (a), (b) and (c) are scattering patterns before shear and at intermediate and late times 
after the onset of shear, respectively, while, for completeness, (d) details the side view of the structure function corresponding to 
(c). The initial condition for shearing was a gyroid on a 128 3 lattice at 15 000 time steps of self-assembling. Time steps after the 
start of shear for these snapshots are indicated below each. Darkness in the greyscale grows with the scattering intensity — filled 
isocurves correspond to values S = 1, 80, 200, 700. The spikes are shear-dependent features; see Fig. 8 and text for discussion. 
All quantities reported are in lattice units, and N = iV x = N z . 
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FIG. 8: Structure function of the 'synthetic gyroid', as calculated using Eq. (7) on the field G(x), the expression of the latter 
being Eq. (9). Parameter 5 max is the maximum value of the dephase <5(x) = (a; — iV x /2)<5 max , which serves to mimic a uniform 
strain across the structure. The case <5 max = gives an approximation to the Schoen G (or 'ideal gyroid') structure. Darkness 
in the greyscale grows with the scattering intensity, and the filled isocurves shown correspond to S = 1, 80, 200, 700. For k z / 0, 
the strain shifts the pattern leftwards and smear the peaks, while leaving the k z = peaks intact. The smearing not being in 
direct relation to the strain — panel (b) shows more smearing than panel (c) — suggests a similar behaviour for the spikes shown 
in Fig. 7(c). The 'cardioid' shape reported in Fig. 7 originates from the facts that the structure undergoes a structural transition 
(weaking and/or relocation of some of its k z / 0, fc x / peaks) whilst being sheared with a velocity profile of positive slope (cf. 
Fig. 2), which orients the 'atria' leftwards. All quantities reported are in lattice units, and N = N x = N z . 
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FIG. 9: Temporal evolution of the average shear stress re- 
sponse of a lamellar mesophase at a shear speed of U — 0.10, 
for different initial amphiphile densities w°^ s . The solid, 
dashed, and dash-dotted curves correspond, respectively, to 
n (o)s _ o gg^ 0.95. The average is computed over the three- 
dimensional domain [8, iV x - 8] x [1, N y ] x [16, N z - 16], 
where jV x = N y = N z = 128 and error bars are not included 
since they are negligible. All quantities reported are in lattice 
units. 
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(a) n<°) s = (b) n<°) s = 0.80 (c) n<°) s = 0.95 



FIG. 10: Slabs < y < 8 of the order parameter <p for sheared lamellar mesophases corresponding to increasing amphiphile 
density, n' -^, as indicated below each relevant panel, at time step At = 8 000 after the onset of shear and for shear velocity 
U = 0.10. The coordinate system is the same as that in Fig. 5. In panel (a), the regions opaque to incident (volume rendering) 
light are those for which <p > 0.18, where < 0.36 across the system. In panel (b), the opaque regions are those for which 
4> > 0.22, where \<p\ < 0.45 across the system. In panel (c), the opaque regions are those for which <j> > 0.24, where \<j>\ < 0.48 
across the system. It is worth noting that the surfactantless case, (a), exhibits a curved interface. The amphiphilic cases, (b) 
and (c), display the formation of irregularities in the interface and nodal planes, as a result of the inter-amphiphile interaction. 
All configurations have translational symmetry along the y-axis. All quantities reported are in lattice units. 



